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ABSTRACT 

The Kozai mechanism often destabihses high incUnation orbits. It couples changes in the 
eccentricity and incHnation, and drives high inclination, circular orbits to low inclination, 
eccentric orbits. I n a recent study of the dynamics of planetesimals in the quadruple star 
system HD98800 (IVerrier & Evansll2008l) . there were significant numbers of stable particles 
in circumbinary polar orbits about the inner binary pair which are apparently able to evade the 
Kozai instability. 

Here, we isolate this feature and investigate the dynamics through numerical and analyt- 
ical models. The results show that the Kozai mechanism of the outer star is disrupted by a 
nodal libration induced by the inner binary pair on a shorter timescale. By empirically mod- 
elling the period of the libration, a criteria for determining the high inclination stability limits 
in general triple systems is derived. The nodal libration feature is interesting and, although ef- 
fecting inclination and node only, shows many parallels to the Kozai mechanism. This raises 
the possibility that high inclination planets and asteroids may be able to survive in multistellar 
systems. 
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1 INTRODUCTION 

The Kozai instability is a well-known destabilizing mechanism of 
high inclination satellites. Arnold (1990), in his book Huygens & 
Barrow, Newton & Hooke, reports a neat illustration due to Lidov 
(1963). Lidov discovered that if the orbit of the Moon is turned 
through 90°, its eccentricity increases so rapidly under the action 
of the tidal forces of the Sun that it collides with the Earth in 
four years! The effect was described and analysed by Kozai (1962, 
1980) in his studies of the survival of high inclination comets and 
asteroids. 

The Kozai mechanism is a secular effect whereby, above a 
critical inclination icrit of around 40° (depending on the system), 
the inclination i is coupled to eccentricity e though the Kozai con- 
stant (see e.g., [Thomas &Morbidellilll996h 



(1) 



and is driven to vary in cycles between its initial value and the crit- 
ical value. The argument of pericentre also librates about ± 90° 
( lKozailll962l : iThomas & Morbi delli' 1996; Takeda et al. 2008). In 
other words, nearly circular, high inclination orbits are driven to 
high eccentricities in exchange for lower inclination, so that in the 
Solar system this produces the population of Sun-grazing comets 
(e.g., Stagg & Bailey 1989). 

At first glance, this suggests that planets or asteroids in 
high inclination orbits in multistellar systems cannot survive for 
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long times. A planet in a circumstellar orbit in a binary star 
system is known to be subject to Kozai cycles on a timescale 
dKiseleva. Eggleton & Mikkolj|l998l ; lTakeda et al]|2008l) 
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where Pbin is the period of the binary, Pp that of the planet, ebin is 
the binary's eccentricity, m,4 and ms the stellar masses and rrip the 
planet's mass. Surprisingly, then, a recent study of planete simals in 
the stellar system of HD 98800 jVerrier & EvaiiibOOSh found a 
long-lived, high inclination circumbinary population. 

HD 98800 has four stars, in two close binaries A and B, in 
highly inclined and eccentric orbits and hosts a circumbinary de- 
bris disc around the B binary pair. The system is well approximated 
as a hierarchical triple star system, as the secondary star in the A 
binary is fairly small and in a close orbit. High inclination parti- 
cles are found to be in long-term stable orbits inclined by 55° to 
135° to the inner binary B. This is well above the critical incli- 
nation relative to both stellar orbits, yet the orbits are not under- 
going Kozai cycles. It is known that the Kozai effect can be dis- 
rupted if another mechanism, such as general relativity, tidal forces 
or other interaction s between planets, causes orbital precession on 
a shorter timescale ( Kinoshita & Nakai 'l99lL IWu & Murravll2003l: 
Takeda & Rasio 2005 ; Takeda etal. 2008). However, the simula- 
tions in lVerrier & Evans! ( |2008|) deal with the Newtonian dynamics 
of test particles, so the only additional factor present is the mutual 
gravitational perturbations of the stars. 

This paper explains the origin of the surprising stability. The 
inner binary causes a nodal precession, instead of Kozai cycles. 
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Figure 1. Test particle distributions for particles starting with initial inclinations in the range 55° to 135° from the simulation using orbit II in lVerrier & Evanj 
120081) . Colour indicates initial outer inclination, as seen in the lower centre panel. Upper and Lower Far Left: The final semimajor axis of survivors plotted 
against the final inclination with respect to the inner and outer binary. Note the survivors are restricted to semimajor axes between 5 and 7 au. Upper and 
Lower Centre Left: Initial and final inclinations of the survivors with respect to the inner and outer binary pair. Upper and Lower Centre Right: The surviving 
distributions in the (i cos Q, i sin Q) plane (relative to the orbital plane and pericentre of the inner binary) show two distinct populations. However, no structure 
is appai'ent in the (e cos u), e sin u)) plane. Upper and Lower Fai" Right: The radial cross sections of the survivors for all three wide orbit simulations: orbit I is 
in black, orbit II in red and orbit III in green. 



Table 1. The initial orbital paramete rs for the simulations of system HD 
98800 used in lVerrier & EvansI <2008l) . These are the orbital elements used 
to evolve the system from 1 Myr in the past to the present state, so the 
starting parameters for the inner binary are also different for each wide orbit 
case. The longitudes given are relative to the arbitrary reference direction 
of the simulations. 



Orbital 


Wide orbit A-B 


Inner orbit Ba-Bb 


Parameter 


I 


II 


III 


I 


II 


III 


mi (Mq) 




1.281 






0.699 




ma (Mq) 




1.3 






0.582 




a (au) 


61.9 


67.6 


78.6 


0.983 


0.983 


0.983 


e 


0.3 


0.5 


0.6 


0.640 


0.790 


0.570 


in 


130.7 


144.7 


131.7 


0.0 


0.0 


0.0 


^(°) 


214.1 


323.4 


65.9 


82.0 


11.3 


34.89 


Q(°) 


224.5 


309.2 


66.6 


0.0 


0.0 


0.0 


M (°) 


166.3 


90.8 


16.2 


225.9 


32.6 


259.7 



This stabilises the test particles against any Kozai instability driven 
by the outer star. We first present the evidence in favour of this 
explanation in Section |2l and then provide a more detailed study 
of the nodal precession, together with the dynamics behind it, in 
Section[3] Finally, we summarize our conclusions in Section|4] 



2 NUMERICAL EXPERIMENTS 

In IVerrier & EvansI ( l2008h , numerical studies of the stability of 
planetesimals in the circumbinary disc of HD9800 were carried us- 
ing the MoiRAI cod e. This is a fast sympl ectic code adapted for 
hierarchical systems l l Verrier & Evan s'2007). We have checked that 
our results are not an artefact of the numerical integration method, 



as similar results are obtained with a standard jPress etalJli98^ 
Bulirsch-Stoer integration scheme (see I Verrieill2008h '. 

The stellar parameters used for the simulations are shown in 
Table [T] We investigate three possible configurations, labelled or- 
bits I, II and III (see Verrier & Evans 2008 for more details on 
the rationale for the orbital parameters). High inclination particles 
are unexpectedly found to be in long-term stable orbits inclined by 
55° to 135° to the inner binary B, and restricted to a narrow range 
of semimajor axis. These particles are shown in Figure [T] for the 
case of the orbit II parameters. Using f2 to denote the longitude of 
the ascending node, we plot the particles in the (i cos SI, i sin Q,) 
plot. This reveals distinct dynamical structures - in particular, two 
populations of stable particles. Interestingly, though, there is little 
structure seen in the (e cos u), e sin uj), where uj is the argument of 
pericentre. 

The orbits of all the stable particles librate about 90° inner 
(i.e. relative to the inner B binary's orbital plane) inclination. Their 
nodes also librate about ±90° relative to the inner binary's pericen- 
tre, defining the two populations seen in Figure [T] These are also 
defined by their initial node: the upper population in the {i cos Q., 
i sin Q.) plot starts at 120° and the lower at 240° (relative to the 
simulations reference direction). All particles starting with Q. = Q° 
are not stable. Unlike the choice of initial node, there appears to 
be no difference in the evolution of particles with differing peri- 
centres or times of pericentre passage. Interestingly, the stability of 
the high inclination particles is greater for the simulations with a 
lower eccentricity for the outer star. This is shown in the far right 
panels of Figure [T] in which there are more surviving test particles 
in the case of orbit I as opposed to orbit III. A typical orbit of one 
of the stable high inclination particles is shown in Figure |2l and 
illustrates the libration of inclination and node. There also appears 



© 0000 RAS, MNRAS 000,[T]l7] 



High Inclination Planets in Multistellar Systems 3 





1.0 








io-90.0000' 


>, 


iliil 


150 


mmmmm 




■t 0.6 




1 100 






§ 0-4 














50 






^ 0.2 
0.0 




- 





0.0 0,2 0-4 0,5 0,8 1-0 
t (Myr) 



0,0 0,2 0.4 0-6 0.8 1-0 
t (Myr) 



0,0 0,2 0,4 0.6 0-8 1.0 
t (Myr) 



300 
200 











-S 300 






200 




t,..A:'\:-." 


g 100 
^ 





0.0 0.2 0.4 0.6 0.8 1-0 
I (Myr) 



0.0 0.2 0-4 0-6 0.3 1.0 
I (Myr) 




Figure 2. The orbital evolution, relative to the inner binary orbit, of a parti- 
cle starting near 7 au with an initial inclination of 90° . While the inclination 
variations over the Myr integration length are large, the eccentricity remains 
small and the particle is clearly not undergoing Kozai cycles. The inner stel- 
lar eccentricity and mutual inclination are oveiplotted as solid lines on the 
relevant panels. 
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Figure 3. The typical short term evolution of a pailicle when the outer star 
is removed. Orbital elements are relative to the inner binary's orbit, using 
the star's periastron as the reference direction. 



to be a slight modulation of this variation with the stellar orbital 
evolution, as indicated in the Figure. 

A brief investigation of the variation of stability with the stel- 
lar parameters was achieved by integrating the stable particles for 
a variety of different eccentricities and mutual inclinations of the 
two stellar orbits. This revealed that the outer stellar eccentricity 
has little effect on the stability, while the inner eccentricity is cru- 
cial! As the inner orbit becomes circular, the high inclination sta- 
bility, somewhat surprisingly, decreases dramatically. The stellar 
mutual inclination is only important insofar as it may cause vari- 
ations in the resulting inner eccentricity. This would appear to be 
in contradiction to the results shown in Figure [T] where the stable 
region increases as the outer stars eccentricity decreases. However, 
in these simulations the test particles do not start similarly aligned 
with the inner star in each case (as the original simulations start at 
1 Myr in the past, where the inner star has a different longitude for 
each outer eccentricity case), meaning that in the Orbit I case the 
particles are closer to the libration island, while in the Orbit III case 
they are further away and hence less stable. This effect is discussed 
in more detail in Section lSTl 

The orbital evolution of particles due to each stellar orbit is 



also of interest. Replacing the inner binary as a single star results 
in the stable particles now undergoing Kozai cycles and becoming 
rapidly unstable at very high inclinations. However, removing the 
outer star instead does not have similar results. In this case, the 
particles remain stable and still librate about 90° inclination and 
longitude of the ascending node, as shown in Figure [3] The libra- 
tions have a period of the order of Kyrs which increases with parti- 
cle semimajor axis and an amplitude that decreases with semimajor 
axis. 

This nodal libration could be suppressing the Kozai cycles 
caused by the outer star, and could explain the distribution of the 
stable particles, if the period is shorter than that of the cycles. At 
3 au, tkoz ~ 3 Kyr and at 8 au tkoz ^ 1-5 Kyr, while the nodal 
libration period at 3 au is about 1 Kyr, indeed shorter than tkoz, 
and at 8au is about 2 Kyr, longer than tkoz, supporting the idea. 

In summary, the numerical evidence suggests that the inner bi- 
nary star causes the nodes and inclinations of particles within about 
7 au to librate on a shorter timescale than the Kozai cycles induced 
by the outer star, thereby suppressing the mechanism and stabil- 
ising their orbits. The simulations also indicate the period of the 
librations increase as the inner binary's eccentricity decreases, ex- 
plaining the lack of stability of the particles in the triple system 
with a circular inner binary. This nodal libration is itself a very in- 
teresting phenomenon, as is the lack of Kozai cycles. 



3 THE NODAL LIBRATION MECHANISM 

3.1 The Elliptic Restricted Three-Body Problem 

The stellar system is now reduced to that of the inner binary only. 
The outer binary can be neglected, as we have just demonstrated 
that the nodal libration takes place in its absence. This arrangement 
is now identical to the elliptic restricted three-body problem. The 
circular restricted three-body problem is the limit Cbin = 0. It is 
most commonly met under the guise of the Copenhagen problem, 
in which the two bodies have equal mass. There has been substan- 
tial effort on the classification of the orbits in the Copenhagen prob- 
lem (see e.g., Contopoulos 1967), but the eccentric case has been 
much less well studied. 

Simulations are run for 10 Kyrs, as the short term behaviour is 
of interest. Since the initial node of the test particles appears to be 
an important factor in the particle's orbital evolution, a set of sim- 
ulations were run looking at a much wider range of this parameter 
The underlying stellar parameters are those corresponding to orbit 
II, whilst the initial test particle grid is shown in Table|2] 

The resulting distributions in inclination and the (i cos Q, 
i sin Q) plane are shown in Figure |4] Once again the two libration 
islands are clear and are at about i = 90° and Q, = ±90° relative 
to the pericentre line of the binary star Figure [5] shows these is- 
lands in more detail for a similar simulation that now has particles 
spaced exactly about these libration centre (i.e. the binary's peri- 
centre has been set to lie along the simulation reference direction) 
and from 0° to 180° inclination. Particles in the centre of the is- 
lands remain there, and the initially low inclination particles also 
remain in low inclination orbits, librating about the centre of the 
(i cos fl, i sin fl) plane while the retrograde particles circulate. This 
libration and critical angle is in striking similarity to Kozai mech- 
anism, where similar libration islands are seen, but for u). Here, 
however, the argument of pericentre still rapidly circulates, and no 
structure is apparent in the (e cos ui, e sin lj) plane. 

A particle's proximity to the centre of the libration island is 
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Table 2. The test particle grid used in the simulations shown in Figures |4] Table 3. Fitted parameters for the period law given in eqn @ 
andia 
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Step Size 
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Inc 111131 1 n 


j\[ 


n 
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0.0029 ± 0.0008 


3.33 ± 0.14 
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0.26 


70° 


0.0023 ± 0.0006 
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0.34 


80° 
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3.31 ± 0.14 


0.45 
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Figure 6. The conservation of the angular um component, L^^, , along the 
inner binary's line of apses for the same particle as Figure|2] 



controlled by its initial node and inclination: higher inclinations 
and nodes place the particle further from the centre and increase 
the period of the variation (which makes them less stable against 
the Kozai cycles of the outer star). In addition, particles that do cir- 
culate in both these simulations show larger eccentricity variations 
and in some cases are not stable. 

This explains why in the original HD 98800 simulations only 
those particles with nodes (relative to the simulations arbitrary ref- 
erence direction) near 120° or, to a lesser extent, 240° survived: the 
other nodal value of 0° placed particles in the circulation regime, 
resulting in destabilising increases in eccentricity. For the orbit I 
case, far more particles are stable as the star's longitude means 
some initially lie very close to the libration centres, while orbit III 
particles are less stable as they are further away from it. 



3.2 An Approximate Integral of Motion 

The surfaces of section (Figures |4] and |5} suggest that there is an 
approximate integral of the motion for the circumbinary test parti- 
cles. From the simulation data, we can verify that this is the case. 
The constant is the test particle's component of angular momentum 
along the direction of the line of apses of the binary star's orbit, 
which is 



Lui, = h sin itp sin Qu. 



(3) 



where itp is the particle's inclination relative to the plane of the 
binary's orbit, and Qui, is the particle's node in the same frame 
of reference and relative to the pericentre of the binary's orbit. An 
example of the conservation of this quantity is given in Figure|6]for 
the test particle from Figure[2l The libration seen in the inclinations 



and nodes of the test particles corresponds to a simple constant 
precession of their orbital planes about this axis. 

This suggests that the physical reason for the test particle's 
orbital evolution is a torque on the orbit due to the component of 
the gravitational force from the binary parallel to the axis of the 
line of apses, labelled as the direction x^i, . This offers an expla- 
nation for the increased stability with eccentricity. As the binary's 
orbit becomes more eccentric it approaches a straight line along the 
direction x^,, , resulting in a greater torque and hence faster preces- 
sion. As the binary's eccentricity decreases, the torque due to forces 
along the star's minor axis (i.e. parallel to the j/^^, direction, using 
the previous notation) will increase, which eventually acts to sup- 
press the precession about the 

3.3 Stability Boundaries 

The period and amplitude of the librations vary with initial inclina- 
tion, node and semimajor axis and the shortest period and minimum 
amplitude occur for particles that start close to the libration centre. 
The point at which the period of these variations equals the Kozai 
timescale of the outer star defines the outer stability boundary. This 
can be modelled empirically by fitting the period of the inclination 
variations. 

The nodal librations are best defined for initial nodes of SI = 
90°. To fit the period, sets of simulations are run with a similar grid 
of test particles as shown in Table[2l but with Q and M fixed at 90° 
and 0° respectively and the semimajor axis varied from 3 to 15 au 
now in steps of 0.5 au. The period of the inclination variations are 
shown in Figure [7] as a function of initial semimajor axis for the 
different inclinations. A power law fits each inclination case very 
well. This gives the following empirical law 



iVa" Kyr 



(4) 



where atp is the test particle semimajor axis and the associated best 
fit parameters and uncertainties are given in Table[3] 

The Kozai timescale for HD 98800 is also shown in Figure[7] 
and meets the 80° inclination line at 7.4 au. This matches up very 
well with the outer edge of the high inclination stable test parti- 
cles for HD 9880 (see the left panels of Figure[TJ. Therefore, this 
accurately describes the stability boundary for HD 98800 and can 
be used as a method of finding a general high inclination stability 
limit. It should also be noted that the minimum eccentricity of 0.45 
obtained by the inner binary in the three star case results in longer 
period variations, over plotted on Figure|7] which probably causes 
the slightly diffuse outer edge of the halo seen in Figure[T] 

We now carry out the same calculation for sets of test particles 
in general binary systems. This empirical model can then be used in 
conjunction with the outer star's Kozai timescale to place limits on 
the high inclination stability in a general hierarchical triple system. 
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Figure 4. The inclination variations and surface of section for test particles around the inner binary only. The orbital elements are all relative to the orbit of 
this binary (i.e. the node is from the pericentre of the binary's orbit). Colour indicates initial longitude of ascending node, with black through blue through 
green through yellow through red representing 0° to 360°. The librationcd islands about f2 = ±90° and i = 90° are clear on the left-hand plot of i cos f2 
versus i sin f2. 
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Figure 5. The inclination variations and surface of section for particles around the inner binary of HD 98800, but now spaced evenly about f2 = 90° and for 
inclinations in the range 0° to 180° . Colour now indicates initial inclination, as shown in the left-hand plot. The initially low inclinations remain in the range 
0° to 30° and are seen in the small libration island about the centre of the right-hand plot. 



Simulations sliow that thie libration islands decrease in size 
as the eccentricity decreases, whilst the critical inclination for the 
librations to occur increases. The mass ratio of the binary, /i = 
msb/imBa + niBb), also has an effect, but not on the geometry 
of the libration islands. Instead, as this ratio decreases from 0.5, the 
particles start to escape from the simple libration patterns previ- 
ously seen. To derive our empirical fit, sets of simulations are again 
run with similar test particle grids, but now with only one inclina- 



tion of 85° to find the longest possible period of variations. The 
stellar mass ratio is varied from 0.1 to 0.5 in steps of 0.1 and the 
eccentricity from 0.1 to 0.8 in steps of 0.1. The binary was gener- 
alised to a separation of 1 au, total mass 1 Mq and with all longi- 
tudes equal to zero. The parameters to the power law were found to 
fit 



f(Kyr) 
-Pbin(yr) 



0.001e-^V-°'f— 1 



(3.37±0.06) 



(5) 
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Figure 7. The period of the inchnation variations plotted as a function of semimajor axis, for test particles around the inner binary of HD 98800. The prograde 
inclinations are shown only, as the retrograde cases are almost identical. 50° initial inclination is in green, 60° in blue, 70° in purple and 80° in black and the 
fits to the data are shown as solid lines. The Kozai timescale of the orbit II outer binary star, shown in red, intersects the period variation lines near the outer 
edge of the orbit 11 stable halo. The period of the inclination variations is also shown in orange for an inner binary with eccentricity e = 0.45, the minimum 
value obtained by the stars in the HD 98800 orbit II triple system, and can be seen to be much longer than the higher eccentricity case. 
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Figure 8. The variation of the amplitude and exponent of the fitted power law with mass ratio and eccentricity for the general binary case. In the top two 
panels different mass ratios are indicated by the colours purple (0.1), blue (0.2), green (0.3), orange (0.4) and red (0.5). In the bottom two panels, different 
eccentricities are indicated by the colours black (0.1), purple (0.2), blue (0.3), light blue (0.4), green (0.5), yellow (0.6), orange (0.7) and red (0.8). The fitted 
lines are also shown. 



This is expe cted to scale with semim ajor axis, as seen in the copla- 
nar limits in I Verrier & EvansI j2007h . To derive this empirical for- 
mula, the powers of e and fj, have been fixed at the values shown 
in the fitting. The amplitude has been fitted (excluding the e = 0.1 
case, which seems poorly modelled by this formula, most likely as 
it is approaching the constant inclination regime), but has an asso- 
ciated error of 0.01. Despite this the equation provides a good fit to 
the data and a low value, as shown in Figure[8] 



4 CONCLUSIONS 

This paper presents a foray into an area of dynamics that has not 
been extensively studied, namely small particles around a highly 
eccentric binary star in a hierarchical triple system. At outset, such 
particles might be expected to be disrupted by the Kozai instability. 

Here, we have demonstrated the existence of stable, high in- 
clination circumbinary test particles. They owe their stability to the 



high eccentricity of the inner binary. This is somewhat surprising, 
as it goes against every expectation of planetary dynamics. The in- 
ner binary, instead of inducing Kozai cycles, causes smooth incli- 
nation variations and nodal precession for certain initial longitudes. 
This suppresses Kozai cycles that would otherwise occur due to the 
outer star in the hierarchical triple. 

An analytical theory to explain these variations is desirable. 
However, the stars are very eccentric, their mass ratio is high, and 
the test particle's inclination is also large, leaving no obvious small 
parameters for the standard perturbation expansions of celestial me- 
chanics. The semimajor axis ratios are not very small, but do not 
expla in why existing secular theories (e.g. Kozai 1962; For d et al] 
I2OOOI) are unable to describe the dynamics. The libration islands, 
inclination variations and critical initial inclination are surprisingly 
familiar to, but contrasting with, the dynamics of the Kozai mech- 
anism. The Kozai mechanism occurs if ti; ~ 0, whereas the nodal 
libration appears to be a similar process occurring when f2 ~ 0. 
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There exists an approximate integral of motion - namely, the test 
particle's component of angular momentum along the direction of 
the line of apses of the binary star's orbit. This also is analogous 
to the Kozai mechanism, which has an approximate integral of mo- 
tion. 

The high inclination libration feature may have important con- 
sequences for planetary stability in circumbinary orbits. Many stel- 
lar members of triples have exchanged into the system, so high mu- 
tual planet-star inclinations are very likely. If there are regions of 
stability then the outlook for planetary systems in these environ- 
ments is more promising than previously thought. 
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